#include "AnaDef.h"
#include "XGLUtils.h"

void PrintCheck(const TString tag, const Double_t xold, const Double_t xnew, const Double_t xbak)
{
  if(xnew != xbak && xbak!=-999){
    printf("\n********** PrintCheck() %s NOT consistent! %.0f %.0f delta: %.0f **********\n", tag.Data(), xnew, xbak, xnew-xbak);
    exit(1);
  }
  if(xold<0)
    printf("%-30s%10.0f\n", tag.Data(), xnew);
  else
    printf("%-30s%10.0f  ->%10.0f (%.3f%%)\n", tag.Data(), xold, xnew, xnew/xold*100);
}

void Scan1DHist(const TH1 *h1, const TString tag)
{
  const Double_t ntot = h1->Integral(0,1000);
  printf("\n");
  for(Int_t ix=0; ix<=h1->GetNbinsX()+1; ix++){
    const Double_t cont = h1->GetBinContent(ix);
    if(!cont)
      continue;
      
    printf("%s (%3.0f): %7.0f / %7.0f = %7.3f%%\n", tag.Data(), h1->GetXaxis()->GetBinLowEdge(ix), cont, ntot, cont/ntot*100);
  }
}

void Scan2DHist(const TH2 *h2, const TString tag)
{
  const Double_t ntot = h2->Integral(0,1000,0,1000);
  printf("\n");
  for(Int_t ix=0; ix<=h2->GetNbinsX()+1; ix++){
    for(Int_t iy=0; iy<=h2->GetNbinsY()+1; iy++){
      const Double_t cont = h2->GetBinContent(ix,iy);
      if(!cont)
        continue;
      
      printf("%s (%3.0f, %3.0f): %7.0f / %7.0f = %7.3f%%\n", tag.Data(), h2->GetXaxis()->GetBinLowEdge(ix), h2->GetYaxis()->GetBinLowEdge(iy), cont, ntot, cont/ntot*100);
    }
  }
}

Double_t GetNevent(const TH2I *hh)
{
  const Int_t ny = hh->GetNbinsY();
  if(hh->Integral(0,1000, ny+1,ny+1)){
    printf("hnpair overflow!!!! %d %0.f\n", ny, hh->Integral(0,1000, ny+1,ny+1));
    Scan2DHist(hh,"ntrk, npair");
    exit(1);
  }

  Double_t neve = 0;
  for(Int_t ii=2; ii<=ny; ii++)
    neve += (ii-1)*hh->Integral(0,10000, ii, ii);

  return neve;
}

void CheckStat()
{
  printf("\nChecking Stat...\n");

  TList *ll=0;
  gDirectory->GetObject("xlu_Hist",ll);

  TH1I * htrig=0x0;
  TH1I * hntrdv1=0X0;
  TH1I * hreuse = 0x0;
  TH1I * hnfillevent = 0x0;
  TH1I * hnfilltrk = 0x0;

  const TString hname[]={"a0_htrig", "a1_hntrdv1", "a26_hreuse", "a6_nfillevent", "a7_nfilltrk"};
  TH1I **hhs[]={&htrig, &hntrdv1, &hreuse, &hnfillevent, &hnfilltrk};
  const Int_t nh=sizeof(hhs)/sizeof(TH1I*);
  const Int_t nname = sizeof(hname)/sizeof(TString);
  if(nh!=nname){
    printf("nh!=nname!! %d %d\n", nh, nname);
    exit(1);
  }
  for(Int_t ii=0; ii<nh; ii++){
    if(ll)
      (*hhs[ii])=(TH1I*)ll->FindObject(hname[ii]);
    else 
      gDirectory->GetObject(hname[ii], *hhs[ii]);

    if(!(*hhs[ii])){
      printf("no hist TH1I!! %d %s\n", ii, hname[ii].Data());
      if(ll)
        ll->ls();
      else
        gDirectory->ls();

      exit(1);
    }
  }

  TH1D * hminphi = 0x0;
  const TString h1n[]={"a5_minphi"};
  TH1D **h1s[]={&hminphi};
  const Int_t n1h = sizeof(h1s)/sizeof(TH1D*);
  const Int_t n1nn = sizeof(h1n)/sizeof(TString);
  if(n1h!=n1nn){
    printf("n1h!=n1nn!! %d %d\n", n1h, n1nn);
    exit(1);
  }
  for(Int_t ii=0; ii<n1h; ii++){
    if(ll)
      (*h1s[ii])=(TH1D*)ll->FindObject(h1n[ii]);
    else 
      gDirectory->GetObject(h1n[ii],*h1s[ii]);
    
    if(!(*h1s[ii])){
      printf("no h1!! TH1D %d %s\n", ii, h1n[ii].Data());
      if(ll)
        ll->ls();
      else
        gDirectory->ls();
      exit(1);
    }
  }

  TH2I * hnpair = 0x0;
  const TString h2in[]={"a25_hnpair"};
  TH2I **h2is[]={&hnpair};
  const Int_t n2ih = sizeof(h2is)/sizeof(TH2I*);
  const Int_t n2inn = sizeof(h2in)/sizeof(TString);
  if(n2ih!=n2inn){
    printf("n2ih!=n2inn!! %d %d\n", n2ih, n2inn);
    exit(1);
  }
  for(Int_t ii=0; ii<n2ih; ii++){
    if(ll)
      (*h2is[ii])=(TH2I*)ll->FindObject(h2in[ii]);
    else 
      gDirectory->GetObject(h2in[ii], *h2is[ii]);
    
    if(!(*h2is[ii])){
      printf("no hist!! TH2I %d %s\n", ii,  h2in[ii].Data());
      if(ll)
        ll->ls();
      else
        gDirectory->ls();
      exit(1);
    }
  }

  //========================================================================

  const Double_t nraw  = htrig->Integral(0,10000);
  const Double_t ntrig = htrig->Integral(2,10000);

  const Double_t prentrd = hntrdv1->Integral(0,10000);
  const Double_t posntrd = hntrdv1->Integral(2,10000);

  const Double_t prenpair = hnpair->Integral(0,10000, 0, 10000);
  const Double_t posnpair = GetNevent(hnpair);
  Scan2DHist(hnpair,"ntrk, npair");

  const Double_t prereuse = hreuse->Integral(0,10000);
  const Double_t posreuse = hreuse->Integral(1,1);

  const Double_t preminphi = hminphi->Integral(0,100000);
  const Double_t posminphi = hminphi->Integral(11,10000);

  const Double_t prenfillevent = hnfillevent->Integral(0,10000);
  const Double_t posnfillevent = hnfillevent->Integral(2,10000);

  const Int_t fillbin1= 12;
  if(hnfilltrk->GetXaxis()->GetBinLowEdge(fillbin1)!=1){
    printf("checkstat fail wrong fillbin1 %d %f\n", fillbin1, hnfilltrk->GetXaxis()->GetBinLowEdge(fillbin1));
    exit(1);
  }

  const Double_t prenfilltrk = hnfilltrk->Integral(0, 1000000);
  const Double_t posnfilltrk = hnfilltrk->Integral(fillbin1, 1000000);
  Scan1DHist(hnfilltrk,"filltrk");
  
  printf("\n");

  PrintCheck("Triggered", nraw, ntrig, prentrd);
  PrintCheck("Ntrdv1", prentrd, posntrd, prenpair);
  PrintCheck("npair > 0", prenpair, posnpair, prereuse);
  PrintCheck("!reuse", prereuse, posreuse, preminphi);
  PrintCheck("minphi > 10", preminphi, posminphi, prenfillevent);
  PrintCheck("xcheck event prefilled", -999, prenfillevent, prenfilltrk/2);
  PrintCheck("xcheck event    filled", prenfillevent, posnfillevent, -999);
  PrintCheck("xcheck track    filled", prenfilltrk, posnfilltrk, -999);

  //============================== soly Drawing...
  TCanvas *c1 = new TCanvas;
  hnpair->SetXTitle("ntrk");
  hnpair->SetYTitle("npair");
  hnpair->SetMarkerSize(2);
  hnpair->GetYaxis()->SetRange(0,10);
  hnpair->Draw("colz");
  gPad->SetLogz();
  XGLUtils::PrintOut("npair", c1, XGLUtils::KPNG);
  delete c1;
}

int main(int argc, char* argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc<1)
    return 1;

  TFile::Open(argv[1]);
  CheckStat();
  return 0;
}
